Single tensionless transition in the Laplacian roughening model 
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We report large scale Monte Carlo simulations of the equilibrium discrete Laplacian roughen- 
ing (dLr) model, originally introduced as the simplest one accommodating the hexatic phase in 
two-dimensional melting. The dLr model is also relevant to surface roughening in Molecular Beam 
Epitaxy (MBE) . We find a single phase transition, possibly of the Kosterlitz-Thouless type, between 
a flat low-temperature phase and a rough, tensionless, high-temperature phase. Thus, earlier con- 
clusions on the order of the phase transition and on the existence of an hexatic phase are seen as 
due to finite size effects, the phase diagram of the dLr model being similar to that of a continuum 
analog previously formulated in the context of surface growth by MBE. 

PACS numbers: 64.60.Cn, 68.35.Ct, 68.35.Rh 



Two-dimensional (2D) melting has played a driving 
role in Statistical Physics for more than two decades. 
Efforts made at clarifying its nature [l| have aided to 
understand systems in which topological defects are rel- 
evant, from the equilibrium fluctuations of metallic sur- 
faces 2] to superfluidity and superconductivity in thin 
films, and phase transitions in liquid crystals Q. One of 
the most intriguing related notions is the hexatic phase, 
between a solid at low temperature (T) and an isotropic 
fluid at high T, transitions between phases being of the 
Kosterlitz-Thouless (KT) type. Such is the Kosterlitz- 
Thoulcss-Halperin-Nelson- Young (KTHNY) mechanism 
for 2D melting Q. Although controversial for some time, 
the hexatic phase has indeed been found in atomistic 
model systems Q and in experiments 

A successful approach to systems with defect-mediated 
phase transitions as the above has been the use of duality 
to formulate equivalent height models. E.g., the discrete 
Gaussian (dG) model [Eq. Q below for bending rigid- 
ity parameter k = 0] is dual of the 2D Coulomb gas, 
and the roughening transition in the former corresponds 
to the well-known KT phase transition of the latter, 
driven by the unbinding of vortex-antivortex pairs. With 
a similar philosophy, the discrete Laplacian roughening 
(dLr) model was introduced by Nelson 7] to describe 2D 
melting. Its Hamiltonian is 

^=^E( CT [ V ^W] 2 + «[V 2 Mr)] 2 ) , (1) 

r 

where r denotes position on a 2D lattice of lateral size 
L, V<j is discrete gradient, and h e Z. The original dLr 
model is obtained by setting to zero the surface tension 
parameter a. Note, the dLr model is a discrete version of 
the linear approximation to Helfrich's energy functional 
for 2D membranes and provides a simplified descrip- 
tion of fluctuating tensionless surfaces, such as biological 
membranes @ or, e.g., such as those grown under condi- 
tions typical in Molecular Beam Epitaxy (MBE) . 



For the dLr model, the KTHNY mechanism would im- 
ply [l| an intermediate hexatic phase in which the surface 
disorders in heights, but not in slopes (quasi-long range 
orientational order). For low T, the surface would be in a 
flat phase, dual of the isotropic fluid in melting, while for 
high T the surface would disorder in heights and slopes, 
providing the dual of the solid phase. In terms of the sur- 
face structure factor S(q) = (h(q)h(—q)) 01 1 the rough 
high T phase implies power law behavior as S(q) ~ q~ 4 , 
changing to S(q) ~ q~ 2 in the hexatic phase jlj, and to 
existence of a finite correlation length in the flat low T 
phase. Equivalently, for the stationary height-difference 
C(r) and slope-difference Cd{r) correlations these 
behaviors amount to: (i) rough phase C(r) ~ r 2 logr, 
Cd{r) ~ logr; (n) hexatic phase C(r) ~ logr, Ca{r) ~ 1; 
(Hi) flat phase C(r) ~ 1, Cd(r) ~ 1. Results supporting 
this picture were obtained on small (L < 32) square and 
triangular lattices 0|. However, conflicting evidence for 
L < 64 was presented that the model had a single first 
order transition, see [l4| and references therein. The dis- 
crepancy has remained unsolved, in spite of recent ana- 
lytical studies [I^i elucidation of the phase diagram be- 
ing important to the diverse contexts mentioned above. 

Here, we provide new Monte Carlo (MC) simulations of 
the dLr model on the square and triangular lattices. Our 
results for sizes up to 512 x 512, much larger than those 
previously studied allow us to see previous works 

as inconclusive due to finite size effects. The model has a 
single continuous transition, possibly of the KT type, be- 
tween the flat and the rough phases, there being no sign 
of an hexatic phase to within our numerical resolution 
in T . Notably, this provides an instance of a roughen- 
ing transition in which the rough phase corresponds to a 
free tensionless surface, rather than a free surface with 
tension, as in the dG model. Moreover, the phase dia- 
gram of the dLr model is seen to resemble closely that 
of a continuum model proposed [9|, Il5| for MBE growth, 
suggesting that both models are in the same universality 
class, much like the relationship between the dG and the 
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FIG. 1: Surface morphologies for three sample temperatures 
around T c on the L = 128 square lattice for Neumann (zero 
derivative) boundary conditions. Inset: Lateral cut of the 
surface for T = 1.7. All units are arbitrary. 




FIG. 2: Effective exponent r in the small wave- vector behav- 
ior S(k) ~ k' r (for k = 2sin(q/2) < k*(L) « 3tt/£), as a 
function of T, for L = 32 (O), 64 (□), and 128 (o). Inset: 
surface structure factor S(k) on the L = 128 square lattice 
vs k, for T = 1.62 (bottom) up to T = 1.69 (top). Dashed 
reference lines have slopes —2 (bottom) and —4 (top). All 
other lines are guides to the eye. All units are arbitrary. 



continuum sine-Gordon models 

For our MC simulations we follow the same procedure 
as in Ha. fluctuations being treated by the histogram 
method |l7j |. further validated through additional simu- 
lations on different points of the extrapolated intervals. 
Thermalization has been checked by monitoring the be- 
havior of non-local observables like the specific heat and 
the structure factor at the smallest wave-vector on our 
finite lattices, S(q = 2%/L), as functions of MC time. 
Note that the dLr model has a richer ground state struc- 
ture than the dG model, Hamiltonian (JIJ with a = 
being minimized not only by configurations with uni- 
form heights, but also by configurations with uniform 
slopes (and by more complex morphologies, see below), 
see e.g. the surface morphology made up of patches with 
various constant slopes shown in Fig. 01 for high T . In 
simulations, this requires large enough system sizes and 
appropriate boundary conditions so that the full minima 
structure can be significatively probed. In particular, for 
small sizes and periodic boundary conditions the system 
is effectively constrained to fluctuating around a single 
energy minimum (the morphology with zero slope), in- 
ducing apparent hystheretic behavior associated with a 
first order transition [T^j . In our simulations, we have 
employed both periodic and free (Neumann) boundary 
conditions, and we have made sure that results provided 
are (qualitatively) independent of these. 

As done for the dG model in [? ], we study the phase 
transition through the behavior of the structure factor 
S(q) for different temperatures. In order to test the 
KTNHY mechanism, we have studied the behavior of 
S(q) as a function of T and L, by fitting the small wave- 
vector part of S(q) to S(q) ~ q~ r . As seen in Fig. Elfor 
the square lattice (for the sake of clarity, we omit plots 
for the triangular lattice, in which completely analogous 



results are obtained), there is no evidence of a finite tem- 
perature interval within which r ~ 2, that would be the 
signature of the hexatic phase. Rather, we find a grad- 
ual change from the flat phase behavior (r ~ 0) to the 
rough phase one (r ~ 4). This change becomes more 
abrupt when the system size is increased, so that only 
the flat and the rough phases remain well-defined in the 
thermodynamic limit. These results may thus explain 
the apparent observation of an hexatic phase in [l3l | for 
small L values, where no systematic finite size effects 
were assessed. By defining the critical temperature T c as 
the value at which curves for different system sizes cross 
[l7|. we estimate T c = 1.65(1) for the square lattice and 
T c = 1.90(2) for the triangular lattice. 

Further evidence on the existence of a single phase 
transition is provided by the behavior the specific heat 
c(T,L) = ((H 2 dhT ) - (HdLr) 2 ) / (T 2 L 2 ) as a function of 
temperature. Fig. ^ shows c(T) on the square and tri- 
angular lattices for the largest system sizes in our simu- 
lations. Within our statistics, a single peak at T = T* 
can be detected, rather than two as would be expected 
within the KTHNY scenario. The height and position 
of the peak are functions of lattice size L. Fig. [21 (inset) 
provides the results of finite size analysis on the specific 
heat curves, in which the maximum value c max (L) ob- 
tained for each lattice size is plotted as a function of L. 
Remarkably, although for lattice sizes L < 70 the spe- 
cific heat grows approximately as c max ^ L — compatible 
with claims on the apparent weakly first order character 
of the transition for L < 64 [Tj| — > for larger L values 
the increase of c max (£) slows down. For our largest sim- 
ulated systems, the best fit is logarithmic c max ~ logi. 
Actually, for the 2D XY model the specific heat at the 
transition temperature is known [Isj to first grow loga- 
rithmically with system size and then saturate for large 
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FIG. 3: c(T) vs T on the square (main panel) and triangular 
(inset) lattices for L = 16 (*), 32 (A), 64 (o), 128 (O), 256 
(□), and 512 (v)- Bars are statistical errors and thin lines are 
guides to the eye. For the larger L values in each case, thick 
solid lines show the c(T) curve extrapolated by the histogram 
method |l7jl . All units are arbitrary. 

enough values of L, suggesting our result might reflect fi- 
nite size effects. Indeed, saturation is expected provided 
that the correlation length at T* is smaller than L and 
thus a horizontal plateau would occur at low k for S{k), 
namely r — as defined in Fig. [2J The steady decrease 
of r at T* for increasing L indicates that such a condi- 
tion has not being reached. Persistence of logarithmic 
behavior in the L — > oo limit would rather suggest that 
the phase transition is in e.g. the 2D Ising class [l^- In 
order to explore this possibility, in Fig. |2 we study the 
dependence of the specific heat jump position T*(L) with 
lateral size L. In a continuous transition, T*(L) scales 
as ,19] T*(L) - T* ~ lr y l v {\ + gL,-"), where g is a nu- 
merical constant and u> is an exponent that accounts for 
corrections to scaling, and is in the range 7/4 < to < 2 for 
the 2D Ising class [2(J. The best multiparameter fit to 
such scaling form yields v = 1.54(47), and v — 0.94(16) 
on the square and the triangular lattices, together with 
uj = 1.6(3), and 2.2(2.0) respectively, to be compared 
with v = 1 for the 2D Ising class |2JJ . Although these 
results might seem compatible with 2D Ising universal- 
ity for the present transition, we believe our numerical 
evidence favors more strongly a different interpretation. 
Thus, in marked contrast with 2D Ising and as shown by 
Fig. |3 the transition in the dLr model is from a phase 
with finite correlation length to a continuous line of fixed 
points [in the Renormalization Group (RG) sense], char- 
acterized by an infinite value of the correlation length, as 
occurs in a KT transition [l^. In order to corroborate 
the latter interpretation, we can try a phenomenological 
KT-type form for T*(L), namely 

T*(L)=T* + 5 , (2) 

(logi + 6) 2 ' 

for constant a and b. As seen in Fig. [21 this fit is in very 



FIG. 4: Transition temperature T* as obtained from Fig.Q as 
a function of L for the square (□) and triangular (A) lattices. 
For each case, the dashed line is a power-law fit T* (L) — T* ~ 
and the solid line is a fit to Eq. (|5J. L = 512 is not 
employed for the fit due to low statistics. Inset: c max (L) vs 
L on the square (□) and triangular (A) lattices. Lines are 
logarithmic fits to the data, shown for reference. All bars 
represent statistical errors and all units are arbitrary. 



good agreement with the numerical data for large sizes. 
We must caution the reader on the well-known feature 
of the KT transition, that the peak of the specific heat 
does not occur at the critical temperature but, rather, 
at a temperature preceding T c 0, [^. Although the 
size of this offset can be model-dependent, Fig. [21 indeed 
provides estimates, T* = 1.63(1) on the square lattice, 
and T* = 1.85(1) on the triangular lattice, that are below 
the corresponding T c values, and are still inside the low T 
behavior for the spatial correlation functions, see Fig. [3] 
Thus, the inexistence of an intermediate phase and the 
fact that the spatial correlations and the specific heat 
change behavior at different values of T can be hardly 
reconciled with a single transition of the Ising class. 

Our results seem to replace the KTHNY scenario for 
the dLr model by a single, tensionless, KT-type phase 
transition. The absence of the hexatic phase may seem 
surprising when contrasted with the often accepted ar- 
gument that, for increasing T, the surface should first 
disorder in heights and, then, in slopes. However, this is 
a sufficient condition for surface roughening, but it is not 
necessary. For instance, in the dG model slopes are not 
disordered at any temperature. It is also conceivable, as 
is our belief, that heights and slopes disorder at the same 
temperature in the dLr model. This remarkable result is 
also against the expectation that discreteness in surface 
heights renormalizes the surface tension cr, as it indeed 
does in the dG model 16]. For the dLr model, gener- 
ation of a non-zero a would imply that the asymptotic 
properties of the high T phase coincide with those of the 
dG model 11]. In order to explore this possibility, var- 
ious analytical approaches [15J have been applied to the 
following continuum analog of the dLr model, introduced 
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in the context of growth by MBE 9] : 

— = - K W i h sin + v/2fc B TC, 3 

at aj_ \ aj_ / 

where £ is a delta-correlated Gaussian white noise, and 
a±, V, are parameters. Although a dynamical RG study 
for Q does predict the generation of a non-zero surface 
tension, numerical simulations of this Langevin equation 
|i| give results in complete qualitative agreement with 
those of the dLr model presented here. The discrepancy 
between the RG arguments and the numerical results for 
both the discrete and continuum models might be due 
to inaccuracies in the treatment of model symmetries in 
the RG studies. Namely, the dLr model can be writ- 
ten as a model for the surface slopes m = V^/i, i.e., 
Ti-dhi — § J2ri^d ■ m(r)] 2 , with the implicit restriction 
that Vd x m = 0. Thus, the dLr model has larger symme- 
tries than the dG model, the Hamiltonian being invariant 
under arbitrary global shifts in the heights, as in the lat- 
ter, but also in the slopes. Thus the ground state degen- 
eracy here is much larger, minima occurring for all height 
configurations for which Vd • m = 0. However, standard 
perturbative RG analyses 0, |2^| are oblivious to such 
an added complexity in the ground state structure of the 
model. Perhaps in a related fashion, the zero-vorticity 
constraint for the slope field may be playing a dynamical 
role in the unbinding of surface defects for T = T c in the 



2D melting transition described by the dLr model. 

Summarizing, we have found that the dLr model fea- 
tures a single continuous phase transition. Although sizes 
of our simulations are confined to a regime in which the 
specific heat still grows logarithmically with L, rather 
than saturating as in proper KT scaling, the combined 
information from the spatial correlations and the specific 
heat are consistent with a KT transition. This behav- 
ior is remarkably similar to that of the continuum model 
@, including the tensionless nature of the high T phase. 
Progress in the analytical description of these phenomena 
might improve our understanding of non-perturbative ef- 
fects in defect-mediated transitions, and of dynamical 
effects of geometrical constraints (such as the curl-free 
condition above) in equilibrium systems. 
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